Feature engineering for electricity load forecasting

Feature engineering for electricity load forecasting#

The purpose of this notebook is to demonstrate how to use skrub and polars to perform feature engineering for electricity load forecasting.

We will build a set of features from different sources:

  • Historical weather data for 10 medium to large urban areas in France;

  • Holidays and calendar features for France;

  • Historical electricity load data for the whole of France.

All these data sources cover a time range from March 23, 2021 to May 31, 2025.

Since our maximum forecasting horizon is 24 hours, we consider that the future weather data is known at a chosen prediction time. Similarly, the holidays and calendar features are known at prediction time for any point in the future.

Therefore, features derived from the weather and calendar data can be used to engineer “future covariates”. Since the load data is our prediction target, we will can also use it to engineer “past covariates” such as lagged features and rolling aggregations.

Environment setup#

We need to install some extra dependencies for this notebook if needed (when running jupyterlite). We need the development version of skrub to be able to use the skrub expressions.

%pip install -q https://pypi.anaconda.org/ogrisel/simple/polars/1.24.0/polars-1.24.0-cp39-abi3-emscripten_3_1_58_wasm32.whl
%pip install -q altair holidays https://pypi.anaconda.org/ogrisel/simple/skrub/0.6.dev0/skrub-0.6.dev0-py3-none-any.whl

ERROR: polars-1.24.0-cp39-abi3-emscripten_3_1_58_wasm32.whl is not a supported wheel on this platform.


Note: you may need to restart the kernel to use updated packages.

Note: you may need to restart the kernel to use updated packages.

# The following 3 imports are only needed to workaround some limitations
# when using polars in a pyodide/jupyterlite notebook.
import tzdata  # noqa: F401
import pandas as pd
from pyarrow.parquet import read_table

import polars as pl
import skrub
from pathlib import Path
import holidays
import warnings

# Ignore warnings from pkg_resources triggered by Python 3.13's multiprocessing.
warnings.filterwarnings("ignore", category=UserWarning, module="pkg_resources")

Time range#

Let’s define a hourly time range from March 23, 2021 to May 31, 2025 that will be used to join the electricity load data and the weather data. The time range is in UTC timezone to avoid any ambiguity when joining with the weather data that is also in UTC.

We wrap the polars dataframe in a skrub variable to benefit from the built-in TableReport display in the notebook. Using the skrub expression system will also be useful later.

time_range_start = pl.datetime(2021, 3, 23, hour=0, time_zone="UTC")
time_range_end = pl.datetime(2025, 5, 31, hour=23, time_zone="UTC")
time = skrub.var(
    "time",
    pl.DataFrame().with_columns(
        pl.datetime_range(
            start=time_range_start,
            end=time_range_end,
            time_zone="UTC",
            interval="1h",
        ).alias("time"),
    ),
)
time
<Var 'time'>
Show graph Var 'time'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook").

To avoid network issues when running this notebook, the necessary data files have already been downloaded and saved in the datasets folder. See the README.md file for instructions to download the data manually if you want to re-run this notebook with more recent data.

data_source_folder = Path("../datasets")
for data_file in sorted(data_source_folder.iterdir()):
    print(data_file)

../datasets/README.md
../datasets/Total Load - Day Ahead _ Actual_202101010000-202201010000.csv
../datasets/Total Load - Day Ahead _ Actual_202201010000-202301010000.csv
../datasets/Total Load - Day Ahead _ Actual_202301010000-202401010000.csv
../datasets/Total Load - Day Ahead _ Actual_202401010000-202501010000.csv
../datasets/Total Load - Day Ahead _ Actual_202501010000-202601010000.csv
../datasets/weather_bayonne.parquet
../datasets/weather_brest.parquet
../datasets/weather_lille.parquet
../datasets/weather_limoges.parquet
../datasets/weather_lyon.parquet
../datasets/weather_marseille.parquet
../datasets/weather_nantes.parquet
../datasets/weather_paris.parquet
../datasets/weather_strasbourg.parquet
../datasets/weather_toulouse.parquet

List of 10 medium to large urban areas to approximately cover most regions in France with a slight focus on most populated regions that are likely to drive electricity demand.

city_names = [
    "paris",
    "lyon",
    "marseille",
    "toulouse",
    "lille",
    "limoges",
    "nantes",
    "strasbourg",
    "brest",
    "bayonne",
]

all_city_weather_raw = {}
for city_name in city_names:
    # all_city_weather_raw[city_name] = skrub.var(
    # f"{city_name}_weather_raw",
    all_city_weather_raw[city_name] = (
        pl.from_arrow(read_table(f"../datasets/weather_{city_name}.parquet"))
    ).with_columns(
        [
            pl.col("time").dt.cast_time_unit(
                "us"
            ),  # Ensure time column has the same type
        ]
    )

all_city_weather_raw["brest"]
shape: (38_688, 7)
timetemperature_2mprecipitationwind_speed_10mcloud_coversoil_moisture_1_to_3cmrelative_humidity_2m
datetime[μs, UTC]f32f32f32f32f32f32
2021-01-01 00:00:00 UTCnullnullnullnullnullnull
2021-01-01 01:00:00 UTCnullnullnullnullnullnull
2021-01-01 02:00:00 UTCnullnullnullnullnullnull
2021-01-01 03:00:00 UTCnullnullnullnullnullnull
2021-01-01 04:00:00 UTCnullnullnullnullnullnull
…………………
2025-05-31 19:00:00 UTC17.51750.012.06945.00.16873.0
2025-05-31 20:00:00 UTC16.26750.09.11447199.00.16877.0
2025-05-31 21:00:00 UTC15.51750.07.55999993.00.16984.0
2025-05-31 22:00:00 UTC15.56750.09.0100.00.1782.0
2025-05-31 23:00:00 UTC15.56750.05.506941100.00.17181.0
all_city_weather_raw["brest"].drop_nulls(subset=["temperature_2m"])
shape: (36_744, 7)
timetemperature_2mprecipitationwind_speed_10mcloud_coversoil_moisture_1_to_3cmrelative_humidity_2m
datetime[μs, UTC]f32f32f32f32f32f32
2021-03-23 00:00:00 UTC4.628null10.086427nullnull94.0
2021-03-23 01:00:00 UTC5.0280.011.1832016.0null95.0
2021-03-23 02:00:00 UTC5.0780.010.9667136.0null94.0
2021-03-23 03:00:00 UTC4.6280.010.4647975.0null93.0
2021-03-23 04:00:00 UTC4.4280.010.4647975.0null92.0
…………………
2025-05-31 19:00:00 UTC17.51750.012.06945.00.16873.0
2025-05-31 20:00:00 UTC16.26750.09.11447199.00.16877.0
2025-05-31 21:00:00 UTC15.51750.07.55999993.00.16984.0
2025-05-31 22:00:00 UTC15.56750.09.0100.00.1782.0
2025-05-31 23:00:00 UTC15.56750.05.506941100.00.17181.0
all_city_weather = time.skb.eval()
for city_name, city_weather_raw in all_city_weather_raw.items():
    all_city_weather = all_city_weather.join(
        city_weather_raw.rename(lambda x: x if x == "time" else x + "_" + city_name),
        on="time",
        how="inner",
    )

all_city_weather = skrub.var(
    "all_city_weather",
    all_city_weather,
)
all_city_weather
<Var 'all_city_weather'>
Show graph Var 'all_city_weather'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook").

Calendar and holidays features#

We leverage the holidays package to enrich the time range with some calendar features such as public holidays in France. We also add some features that are useful for time series forecasting such as the day of the week, the day of the year, and the hour of the day.

Note that the holidays package requires us to extract the date for the French timezone.

Similarly for the calendar features: all the time features are extracted from the time in the French timezone.

holidays_fr = holidays.country_holidays("FR", years=range(2019, 2026))

fr_time = pl.col("time").dt.convert_time_zone("Europe/Paris")
calendar = time.with_columns(
    [
        fr_time.dt.date().is_in(holidays_fr.keys()).alias("is_holiday_fr"),
        fr_time.dt.weekday().alias("day_of_week_fr"),
        fr_time.dt.ordinal_day().alias("day_of_year_fr"),
        fr_time.dt.hour().alias("hour_of_day_fr"),
    ],
)
calendar
<CallMethod 'with_columns'>
Show graph Var 'time' CallMethod 'with_columns'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook").

Electricity load data#

Finally we load the electricity load data. This data will both be used as a target variable but also to craft some lagged and window-aggregated features.

load_data_files = [
    data_file
    for data_file in sorted(data_source_folder.iterdir())
    if data_file.name.startswith("Total Load - Day Ahead")
    and data_file.name.endswith(".csv")
]

electricity_raw = skrub.var(
    "electricity_raw",
    pl.concat(
        [
            pl.from_pandas(pd.read_csv(data_file, na_values=["N/A", "-"])).drop(
                ["Day-ahead Total Load Forecast [MW] - BZN|FR"]
            )
            for data_file in load_data_files
        ],
        how="vertical",
    ),
)
electricity_raw
<Var 'electricity_raw'>
Show graph Var 'electricity_raw'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook"). 

electricity = (
    electricity_raw.with_columns(
        [
            pl.col("Time (UTC)")
            .str.split(by=" - ")
            .list.first()
            .str.to_datetime("%d.%m.%Y %H:%M", time_zone="UTC")
            .alias("time"),
        ]
    )
    .drop(["Time (UTC)"])
    .rename({"Actual Total Load [MW] - BZN|FR": "load_mw"})
    .filter(pl.col("time").dt.minute().eq(0))
    .filter(pl.col("time") >= time_range_start)
    .filter(pl.col("time") <= time_range_end)
    .select(["time", "load_mw"])
)
electricity
<CallMethod 'select'>
Show graph Var 'electricity_raw' CallMethod 'with_columns' CallMethod 'drop' CallMethod 'rename' CallMethod 'filter' CallMethod 'filter' CallMethod 'filter' CallMethod 'select'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook"). 

electricity.filter(pl.col("load_mw").is_null())
<CallMethod 'filter'>
Show graph Var 'electricity_raw' CallMethod 'with_columns' CallMethod 'drop' CallMethod 'rename' CallMethod 'filter' CallMethod 'filter' CallMethod 'filter' CallMethod 'select' CallMethod 'filter'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook"). 

electricity.filter(
    (pl.col("time") > pl.datetime(2021, 10, 30, hour=10, time_zone="UTC"))
    & (pl.col("time") < pl.datetime(2021, 10, 31, hour=10, time_zone="UTC"))
).skb.eval().plot.line(x="time:T", y="load_mw:Q")

electricity = electricity.with_columns([pl.col("load_mw").interpolate()])
electricity.filter(
    (pl.col("time") > pl.datetime(2021, 10, 30, hour=10, time_zone="UTC"))
    & (pl.col("time") < pl.datetime(2021, 10, 31, hour=10, time_zone="UTC"))
).skb.eval().plot.line(x="time:T", y="load_mw:Q")

Check that the number of rows matches our expectations based on the number of hours that separate the first and the last dates. We can do that by joining with the time range dataframe and checking that the number of rows stays the same.

assert (
    time.join(electricity, on="time", how="inner").shape[0] == time.shape[0]
).skb.eval()

Lagged features#

We can now create some lagged features from the electricity load data.

We will create 3 hourly lagged features, 1 daily lagged feature, and 1 weekly lagged feature. We will also create a rolling median and inter-quartile feature over the last 24 hours and over the last 7 days.

def iqr(col, *, window_size: int):
    """Inter-quartile range (IQR) of a column."""
    return col.rolling_quantile(0.75, window_size=window_size) - col.rolling_quantile(
        0.25, window_size=window_size
    )


electricity_lagged = electricity.with_columns(
    [pl.col("load_mw").shift(i).alias(f"load_mw_lag_{i}h") for i in range(1, 4)]
    + [
        pl.col("load_mw").shift(24).alias("load_mw_lag_1d"),
        pl.col("load_mw").shift(24 * 7).alias("load_mw_lag_1w"),
        pl.col("load_mw")
        .rolling_median(window_size=24)
        .alias("load_mw_rolling_median_24h"),
        pl.col("load_mw")
        .rolling_median(window_size=24 * 7)
        .alias("load_mw_rolling_median_7d"),
        iqr(pl.col("load_mw"), window_size=24).alias("load_mw_iqr_24h"),
        iqr(pl.col("load_mw"), window_size=24 * 7).alias("load_mw_iqr_7d"),
    ],
)
electricity_lagged
<CallMethod 'with_columns'>
Show graph Var 'electricity_raw' CallMethod 'with_columns' CallMethod 'drop' CallMethod 'rename' CallMethod 'filter' CallMethod 'filter' CallMethod 'filter' CallMethod 'select' CallMethod 'with_columns' CallMethod 'with_columns'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook"). 

import altair


altair.Chart(electricity_lagged.tail(100).skb.eval()).transform_fold(
    [
        "load_mw",
        "load_mw_lag_1h",
        "load_mw_lag_2h",
        "load_mw_lag_3h",
        "load_mw_lag_1d",
        "load_mw_lag_1w",
        "load_mw_rolling_median_24h",
        "load_mw_rolling_median_7d",
        "load_mw_iqr_24h",
        "load_mw_iqr_7d",
    ],
    as_=["key", "load_mw"],
).mark_line(tooltip=True).encode(x="time:T", y="load_mw:Q", color="key:N").interactive()

Investigating outliers in the lagged features#

Let’s use the skrub.TableReport tool to look at the plots of the marginal distribution of the lagged features.

from skrub import TableReport

TableReport(electricity_lagged.skb.eval())

Processing column   1 / 11

Processing column   2 / 11

Processing column   3 / 11

Processing column   4 / 11

Processing column   5 / 11

Processing column   6 / 11

Processing column   7 / 11

Processing column   8 / 11

Processing column   9 / 11

Processing column  10 / 11

Processing column  11 / 11

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook").

Let’s extract the dates where the inter-quartile range of the load is greater than 15,000 MW.

electricity_lagged.filter(pl.col("load_mw_iqr_7d") > 15_000)[
    "time"
].dt.date().unique().sort().to_list().skb.eval()

[datetime.date(2021, 12, 26),
 datetime.date(2021, 12, 27),
 datetime.date(2021, 12, 28),
 datetime.date(2022, 1, 7),
 datetime.date(2022, 1, 8),
 datetime.date(2023, 1, 19),
 datetime.date(2023, 1, 20),
 datetime.date(2023, 1, 21),
 datetime.date(2024, 1, 10),
 datetime.date(2024, 1, 11),
 datetime.date(2024, 1, 12),
 datetime.date(2024, 1, 13)]

We observe 3 date ranges with high inter-quartile range. Let’s plot the electricity load and the lagged features for the first data range along with the weather data for Paris.

altair.Chart(
    electricity_lagged.filter(
        (pl.col("time") > pl.datetime(2021, 12, 1, time_zone="UTC"))
        & (pl.col("time") < pl.datetime(2021, 12, 31, time_zone="UTC"))
    ).skb.eval()
).transform_fold(
    [
        "load_mw",
        "load_mw_iqr_7d",
    ],
).mark_line(
    tooltip=True
).encode(
    x="time:T", y="value:Q", color="key:N"
).interactive()

altair.Chart(
    all_city_weather.filter(
        (pl.col("time") > pl.datetime(2021, 12, 1, time_zone="UTC"))
        & (pl.col("time") < pl.datetime(2021, 12, 31, time_zone="UTC"))
    ).skb.eval()
).transform_fold(
    [f"temperature_2m_{city_name}" for city_name in city_names],
).mark_line(
    tooltip=True
).encode(
    x="time:T", y="value:Q", color="key:N"
).interactive()

Based on the plots above, we can see that the electricity load was high just before the Christmas holidays due to low temperatures. Then the load suddenly dropped because temperatures went higher right at the start of the end-of-year holidays.

So those outliers do not seem to be caused to a data quality issue but rather due to a real change in the electricity load demand. We could conduct similar analysis for the other date ranges with high inter-quartile range but we will skip that for now.

If we had observed significant data quality issues over extended periods of time could have been addressed by removing the corresponding rows from the dataset. However, this would make the lagged and windowing feature engineering challenging to reimplement correctly. A better approach would be to keep a contiguous dataset assign 0 weights to the affected rows when fitting or evaluating the trained models via the use of the sample_weight parameter.

Final dataset#

We now assemble the dataset that will be used to train and evaluate the forecasting models via backtesting.

prediction_time = time = skrub.var(
    "prediction_time",
    pl.DataFrame().with_columns(
        pl.datetime_range(
            start=time_range_start + pl.duration(days=7),
            end=time_range_end - pl.duration(hours=24),
            time_zone="UTC",
            interval="1h",
        ).alias("prediction_time"),
    ),
)
prediction_time
<Var 'prediction_time'>
Show graph Var 'prediction_time'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook"). 

features = (
    (
        prediction_time.join(
            electricity_lagged, left_on="prediction_time", right_on="time"
        )
        .join(all_city_weather, left_on="prediction_time", right_on="time")
        .join(calendar, left_on="prediction_time", right_on="time")
    )
    .drop("prediction_time")
    .skb.mark_as_X()
)
features
<CallMethod 'drop'>
Show graph Var 'prediction_time' CallMethod 'join' Var 'electricity_raw' CallMethod 'with_columns' CallMethod 'drop' CallMethod 'rename' CallMethod 'filter' CallMethod 'filter' CallMethod 'filter' CallMethod 'select' CallMethod 'with_columns' CallMethod 'with_columns' CallMethod 'join' Var 'all_city_weather' CallMethod 'join' Var 'time' CallMethod 'with_columns' X: CallMethod 'drop'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook"). 

horizons = range(1, 25)  # Forecasting horizons from 1 to 24 hours
target_column_name_pattern = "load_mw_horizon_{horizon}h"

targets = prediction_time.join(
    electricity.with_columns(
        [
            pl.col("load_mw")
            .shift(-h)
            .alias(target_column_name_pattern.format(horizon=h))
            for h in horizons
        ]
    ),
    left_on="prediction_time",
    right_on="time",
)

horizon_of_interest = 2
target_column_name = target_column_name_pattern.format(horizon=horizon_of_interest)
predicted_target_column_name = "predicted_" + target_column_name
target = targets[target_column_name].skb.mark_as_y()

from sklearn.ensemble import HistGradientBoostingRegressor


predictions = features.skb.apply(
    HistGradientBoostingRegressor(
        random_state=0,
        learning_rate=skrub.choose_float(
            0.01, 0.9, default=0.1, log=True, name="learning_rate"
        ),
        max_leaf_nodes=skrub.choose_int(
            3, 300, default=30, log=True, name="max_leaf_nodes"
        ),
    ),
    y=target,
)
predictions
<Apply HistGradientBoostingRegressor>
Show graph Var 'prediction_time' CallMethod 'join' CallMethod 'join' Var 'electricity_raw' CallMethod 'with_columns' CallMethod 'drop' CallMethod 'rename' CallMethod 'filter' CallMethod 'filter' CallMethod 'filter' CallMethod 'select' CallMethod 'with_columns' CallMethod 'with_columns' CallMethod 'with_columns' CallMethod 'join' Var 'all_city_weather' CallMethod 'join' Var 'time' CallMethod 'with_columns' X: CallMethod 'drop' Apply HistGradientBoostingRegressor y: GetItem 'load_mw_horizon_2h'

Result:

Please enable javascript

The skrub table reports need javascript to display correctly. If you are displaying a report in a Jupyter notebook and you see this message, you may need to re-execute the cell or to trust the notebook (button on the top right or "File > Trust notebook"). 

altair.Chart(
    pl.concat(
        [
            targets.skb.eval(),
            predictions.rename(
                {target_column_name: predicted_target_column_name}
            ).skb.eval(),
        ],
        how="horizontal",
    ).tail(24 * 7)
).transform_fold(
    [target_column_name, predicted_target_column_name],
).mark_line(
    tooltip=True
).encode(
    x="prediction_time:T", y="value:Q", color="key:N"
).interactive()

from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import make_scorer, mean_absolute_percentage_error, get_scorer

mape_scorer = make_scorer(mean_absolute_percentage_error)

ts_cv_5 = TimeSeriesSplit(n_splits=5, max_train_size=10_000, gap=24)

predictions.skb.cross_validate(
    cv=ts_cv_5,
    scoring={
        "r2": get_scorer("r2"),
        "mape": mape_scorer,
    },
    return_train_score=True,
    verbose=1,
    n_jobs=-1,
).round(3)

[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.

[Parallel(n_jobs=-1)]: Done   5 out of   5 | elapsed:    5.2s finished
fit_time score_time test_r2 train_r2 test_mape train_mape
0 1.088 0.066 0.970 0.997 0.021 0.008
1 1.365 0.065 0.983 0.996 0.020 0.010
2 1.431 0.061 0.972 0.995 0.025 0.010
3 1.389 0.066 0.973 0.994 0.023 0.012
4 1.039 0.039 0.977 0.993 0.023 0.014 
ts_cv_3 = TimeSeriesSplit(n_splits=3, max_train_size=10_000, gap=24)
randomized_search = predictions.skb.get_randomized_search(
    cv=ts_cv_3,
    scoring="r2",
    n_iter=30,
    fitted=True,
    verbose=1,
    n_jobs=-1,
)
randomized_search.results_

Fitting 3 folds for each of 30 candidates, totalling 90 fits
learning_rate max_leaf_nodes mean_test_score
0 0.109154 93 0.980190
1 0.169236 32 0.979379
2 0.142631 199 0.978931
3 0.056970 276 0.978662
4 0.146226 232 0.978143
5 0.167918 224 0.977949
6 0.235284 6 0.975711
7 0.320975 23 0.975268
8 0.035245 42 0.974230
9 0.397895 169 0.973812
10 0.611530 4 0.971484
11 0.104043 5 0.971048
12 0.134697 4 0.970058
13 0.559689 145 0.969951
14 0.730922 11 0.967933
15 0.338535 3 0.966256
16 0.157063 3 0.964471
17 0.675817 36 0.964043
18 0.029972 21 0.963098
19 0.711032 21 0.962959
20 0.036459 9 0.957432
21 0.843862 98 0.955367
22 0.020217 68 0.951587
23 0.019238 132 0.950949
24 0.032555 8 0.948631
25 0.019849 14 0.923808
26 0.015802 15 0.892188
27 0.029450 3 0.888511
28 0.020669 4 0.869532
29 0.010737 33 0.826310 
randomized_search.plot_results()

# nested_cv_results = skrub.cross_validate(
#     environment=predictions.skb.get_data(),
#     pipeline=randomized_search,
#     cv=ts_cv_5,
#     scoring={
#         "r2": get_scorer("r2"),
#         "mape": mape_scorer,
#     },
#     n_jobs=-1,
#     return_pipeline=True,
# ).round(3)
# nested_cv_results

# for outer_cv_idx in range(len(nested_cv_results)):
#     print(
#         nested_cv_results.loc[outer_cv_idx, "pipeline"]
#         .results_.loc[0]
#         .round(3)
#         .to_dict()
#     )

# from joblib import Parallel, delayed

# cv_predictions = []
# for ts_cv_train_idx, ts_cv_test_idx in ts_cv_5.split(prediction_time.skb.eval()):
#     features[ts_cv_train_idx].fit

predictions = features.skb.apply(
    skrub.SelectCols(
        cols=skrub.choose_from(
            [
                skrub.selectors.all(),
                skrub.selectors.filter_names(
                    lambda name: name.startswith("load_mw_"),
                ),
                skrub.selectors.filter_names(
                    lambda name: name.startswith("load_mw_")
                    or name.startswith("temperature_"),
                ),
                skrub.selectors.filter_names(
                    lambda name: name.startswith("load_mw_")
                    or name.startswith("precipitation_"),
                ),
                skrub.selectors.filter_names(
                    lambda name: name.startswith("load_mw_")
                    or name.startswith("wind_speed_"),
                ),
                skrub.selectors.filter_names(
                    lambda name: name.startswith("load_mw_")
                    or name.startswith("cloud_cover_"),
                ),
                # calendar features
                skrub.selectors.filter_names(
                    lambda name: name.startswith("load_mw_") or name.endswith("_fr"),
                ),
            ],
            name="feature_subset",
        )
    )
).skb.apply(
    HistGradientBoostingRegressor(
        random_state=0,
        learning_rate=skrub.choose_float(
            0.01, 0.9, default=0.1, log=True, name="learning_rate"
        ),
    ),
    y=target,
)

randomized_search = predictions.skb.get_randomized_search(
    cv=ts_cv_3,
    scoring="r2",
    n_iter=30,
    fitted=True,
    verbose=1,
    n_jobs=-1,
)
randomized_search.results_

Fitting 3 folds for each of 30 candidates, totalling 90 fits
feature_subset learning_rate mean_test_score
0 all() 0.150188 0.979184
1 all() 0.041626 0.975422
2 all() 0.035388 0.972409
3 filter_names(<lambda>) 0.224432 0.968202
4 filter_names(<lambda>) 0.424321 0.963128
5 all() 0.752968 0.961317
6 filter_names(<lambda>) 0.486231 0.961032
7 all() 0.852650 0.955723
8 filter_names(<lambda>) 0.017859 0.913571
9 filter_names(<lambda>) 0.124354 0.901224
10 filter_names(<lambda>) 0.054564 0.901186
11 filter_names(<lambda>) 0.069022 0.900651
12 filter_names(<lambda>) 0.159816 0.900574
13 filter_names(<lambda>) 0.069133 0.900144
14 filter_names(<lambda>) 0.084425 0.898697
15 filter_names(<lambda>) 0.050684 0.898049
16 filter_names(<lambda>) 0.042047 0.897698
17 filter_names(<lambda>) 0.035197 0.896331
18 filter_names(<lambda>) 0.281298 0.889728
19 filter_names(<lambda>) 0.415160 0.888359
20 filter_names(<lambda>) 0.464677 0.884654
21 filter_names(<lambda>) 0.533526 0.882740
22 filter_names(<lambda>) 0.390743 0.879096
23 filter_names(<lambda>) 0.540765 0.866948
24 all() 0.012603 0.865275
25 filter_names(<lambda>) 0.782276 0.858789
26 filter_names(<lambda>) 0.799087 0.855585
27 filter_names(<lambda>) 0.016800 0.854399
28 filter_names(<lambda>) 0.013759 0.819089
29 filter_names(<lambda>) 0.011550 0.783367 
# nested_cv_results = skrub.cross_validate(
#     environment=predictions.skb.get_data(),
#     pipeline=randomized_search,
#     cv=ts_cv_5,
#     scoring={
#         "r2": get_scorer("r2"),
#         "mape": mape_scorer,
#     },
#     n_jobs=-1,
#     return_pipeline=True,
# ).round(3)
# nested_cv_results

from sklearn.multioutput import MultiOutputRegressor

model = MultiOutputRegressor(
    estimator=HistGradientBoostingRegressor(
        random_state=0,
        learning_rate=skrub.choose_float(
            0.01, 0.9, default=0.1, log=True, name="learning_rate"
        ),
    ),
)

predictions = features.skb.apply(
    model, y=targets.skb.drop(cols=["prediction_time", "load_mw"]).skb.mark_as_y()
).skb.set_name("model")

from sklearn.metrics import r2_score


def multioutput_scorer(regressor, X, y, score_func, score_name):
    y_pred = regressor.predict(X)
    return {
        f"{score_name}_horizon_{h}h": score
        for h, score in enumerate(
            score_func(y, y_pred, multioutput="raw_values"), start=1
        )
    }


def scoring(regressor, X, y):
    return {
        **multioutput_scorer(regressor, X, y, mean_absolute_percentage_error, "mape"),
        **multioutput_scorer(regressor, X, y, r2_score, "r2"),
    }


predictions.skb.cross_validate(
    cv=ts_cv_5,
    scoring=scoring,
    return_train_score=True,
    verbose=1,
    n_jobs=-1,
).round(3)

[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.

[Parallel(n_jobs=-1)]: Done   5 out of   5 | elapsed:   59.5s finished
fit_time score_time test_mape_horizon_1h train_mape_horizon_1h test_mape_horizon_2h train_mape_horizon_2h test_mape_horizon_3h train_mape_horizon_3h test_mape_horizon_4h train_mape_horizon_4h ... test_r2_horizon_20h train_r2_horizon_20h test_r2_horizon_21h train_r2_horizon_21h test_r2_horizon_22h train_r2_horizon_22h test_r2_horizon_23h train_r2_horizon_23h test_r2_horizon_24h train_r2_horizon_24h
0 26.136 2.463 0.013 0.006 0.020 0.008 0.028 0.009 0.036 0.009 ... 0.855 0.994 0.896 0.995 0.873 0.995 0.902 0.995 0.890 0.995
1 32.594 2.350 0.012 0.007 0.020 0.010 0.026 0.011 0.029 0.012 ... 0.943 0.993 0.946 0.994 0.946 0.994 0.947 0.994 0.945 0.994
2 34.048 2.364 0.016 0.007 0.025 0.010 0.030 0.012 0.033 0.012 ... 0.898 0.991 0.918 0.992 0.920 0.993 0.924 0.993 0.922 0.993
3 33.446 2.458 0.015 0.009 0.023 0.012 0.026 0.014 0.027 0.014 ... 0.950 0.991 0.958 0.992 0.958 0.992 0.957 0.992 0.956 0.993
4 24.395 1.388 0.014 0.009 0.023 0.013 0.028 0.014 0.032 0.014 ... 0.934 0.990 0.939 0.993 0.941 0.993 0.940 0.993 0.942 0.993

5 rows × 98 columns

learner = predictions.skb.get_pipeline(fitted=True)

sklearn_learner = learner.find_fitted_estimator("model")
sklearn_learner.estimators_

[HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0),
 HistGradientBoostingRegressor(random_state=0)]

target_column_names = [target_column_name_pattern.format(horizon=h) for h in horizons]
predicted_target_column_names = [
    f"predicted_{target_column_name}" for target_column_name in target_column_names
]

output = pl.concat(
    [
        targets.skb.select(cols=["prediction_time", "load_mw"]).skb.eval(),
        predictions.rename(
            {k: v for k, v in zip(target_column_names, predicted_target_column_names)}
        ).skb.eval(),
    ],
    how="horizontal",
).tail(24 * 14)

import datetime


def plot_horizon(output, prediction_time, past_true_values=5):
    # Find the row that contains the forecasted values
    row = output.filter(pl.col("prediction_time") == prediction_time).rows()[0]
    # Build the forecast dataframe by concatenating the different horizons
    forecast = pl.DataFrame(
        {
            "forecast_time": row[0] + datetime.timedelta(hours=h),
            "forecast_value": predicted_value,
        }
        for h, predicted_value in enumerate(row[2:], start=1)
    )
    # Get the past true values and also the future true values of to the last horizon
    true_values = output.filter(
        (
            pl.col("prediction_time")
            >= prediction_time - datetime.timedelta(hours=past_true_values)
        )
        & (
            pl.col("prediction_time")
            <= prediction_time + datetime.timedelta(hours=len(row) - 2)
        )
    )
    # plot the true and forecasted values
    true_values_chart = (
        altair.Chart(true_values)
        .transform_fold(["load_mw"])
        .mark_line(tooltip=True)
        .encode(x="prediction_time:T", y="load_mw:Q", color="key:N")
    )
    forecast_chart = (
        altair.Chart(forecast)
        .transform_fold(["forecast_value"])
        .mark_line(tooltip=True)
        .encode(x="forecast_time:T", y="forecast_value:Q", color="key:N")
    )
    return (true_values_chart + forecast_chart).interactive()

time_to_plot = datetime.datetime(2025, 5, 25, 0, 0, tzinfo=datetime.timezone.utc)
plot_horizon(output, time_to_plot, past_true_values=24 * 5)
Contents